We need the sf package:
> if (! "sf" %in% rownames(installed.packages())) install.packages("sf")
> library(sf)
Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
The coordinates of the houses of the Hanam cohort are in a zipping shapefile that can be downloaded as so:
> download.file("https://www.dropbox.com/s/yyg6gx0nycpv60b/HH%20point%20Hanam.zip?raw=1", "hanam.zip")
Once downloaded, unzip the file:
> unzip("hanam.zip")
Once unzipped read the shapefile:
> hanam <- sf::st_read("Household_point.shp")
Reading layer `Household_point' from data source `/Users/choisy/Dropbox/aaa/projects/Hanam/Household_point.shp' using driver `ESRI Shapefile'
Simple feature collection with 263 features and 22 fields
geometry type: POINT
dimension: XY
bbox: xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
And clean the disk
> file.remove(grep("^Household", dir(), value = TRUE))
[1] TRUE TRUE TRUE TRUE TRUE TRUE
and
> file.remove("hanam.zip")
[1] TRUE
Here the data, which is of sf class
> hanam
Simple feature collection with 263 features and 22 fields
geometry type: POINT
dimension: XY
bbox: xmin: 105.9174 ymin: 20.48691 xmax: 105.9356 ymax: 20.51665
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
First 10 features:
Header Name Descriptio Type Position tempX X Y Altitude Depth Proximity Temperatur
1 Waypoint H001 14-JAN-09 9:20:51AM User Waypoint N20.49189 E105.92737 N20.49189 20.49189 105.9274 314 m <NA> <NA> <NA>
2 Waypoint H002 14-JAN-09 9:18:10AM User Waypoint N20.49196 E105.92861 N20.49196 20.49196 105.9286 318 m <NA> <NA> <NA>
3 Waypoint H003 <NA> User Waypoint N20.49663 E105.93138 N20.49663 20.49663 105.9314 0 m <NA> <NA> <NA>
4 Waypoint H004 14-JAN-09 9:20:08AM User Waypoint N20.49191 E105.92738 N20.49191 20.49191 105.9274 317 m <NA> <NA> <NA>
5 Waypoint H005 <NA> User Waypoint N20.49529 E105.93026 N20.49529 20.49529 105.9303 6 m <NA> <NA> <NA>
6 Waypoint H006 <NA> User Waypoint N20.49530 E105.93292 N20.49530 20.49530 105.9329 -3 m <NA> <NA> <NA>
7 Waypoint H007 14-JAN-09 2:10:03PM User Waypoint N20.49722 E105.92956 N20.49722 20.49722 105.9296 1 m <NA> <NA> <NA>
8 Waypoint H009 14-JAN-09 3:03:29PM User Waypoint N20.49759 E105.92656 N20.49759 20.49759 105.9266 1 m <NA> <NA> <NA>
9 Waypoint H010 <NA> User Waypoint N20.49672 E105.92964 N20.49672 20.49672 105.9296 -3 m <NA> <NA> <NA>
10 Waypoint H011 <NA> User Waypoint N20.49684 E105.93064 N20.49684 20.49684 105.9306 2 m <NA> <NA> <NA>
Display_Mo Color Symbol Facility City State Country Date_Modif Link Categories geometry
1 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9274 20.49189)
2 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9286 20.49196)
3 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9314 20.49663)
4 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9274 20.49191)
5 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9303 20.49529)
6 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9329 20.4953)
7 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9296 20.49722)
8 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9266 20.49759)
9 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9296 20.49672)
10 Symbol & Name Black Residence <NA> <NA> <NA> <NA> <NA> <NA> <NA> POINT (105.9306 20.49684)
The coordinates of the houses can be extracted from that object and converted into a data frame as so:
> houses <- as.data.frame(hanam[, c("Name", "X", "Y")])[, -4]
which gives:
> head(houses)
Name X Y
1 H001 20.49189 105.9274
2 H002 20.49196 105.9286
3 H003 20.49663 105.9314
4 H004 20.49191 105.9274
5 H005 20.49529 105.9303
6 H006 20.49530 105.9329
We need the OpenStreetMap package:
> if (! "sf" %in% rownames(installed.packages())) install.packages("OpenStreetMap")
> library(OpenStreetMap)
Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
Below is the code you’d need to download the tiles from 3 different types of maps:
> upperleft <- c(20.525681, 105.905383)
> lowerright <- c(20.463343, 105.969019)
> bing <- openmap(upperleft, lowerright, type = "bing", minNumTiles = 20)
> osm <- openmap(upperleft, lowerright, type = "osm", minNumTiles = 20)
> esri <- openmap(upperleft, lowerright, type = "esri", minNumTiles = 20)
But the downloading would be twice as fast if you downlaod the objects directly from here:
> download.file("https://www.dropbox.com/s/26y2pgouodj4rpe/bing.rds?raw=1", "bing_file.rds")
> download.file("https://www.dropbox.com/s/rd34k3ixwthzg9g/esri.rds?raw=1", "esri_file.rds")
> download.file("https://www.dropbox.com/s/827ze3xr0orbab0/osm.rds?raw=1" , "osm_file.rds")
> bing <- readRDS("BING.rds")
> esri <- readRDS("ESRI.rds")
> osm <- readRDS("OSM.rds")
> for(file in paste0(c("bing", "esri", "osm"), "_file.rds")) file.remove(file)
The function that makes the map:
> plot_hh <- function(map, points) {
+ require(sf)
+ plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)))
+ plot(map, add = TRUE)
+ plot(st_geometry(st_transform(points, map$tiles[[1]]$projection)), add = TRUE, col = "red")
+ }
OSM map:
> plot_hh(osm, hanam)
ESRI map:
> plot_hh(esri, hanam)
BING map:
> plot_hh(bing, hanam)